home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / dsp / dspgroup / fir.arc / FIR.FOR < prev   
Encoding:
Text File  |  1987-08-01  |  19.5 KB  |  701 lines
  1. $NOFLOATCALLS
  2. C
  3. C-----------------------------------------------------------------------
  4. C MAIN PROGRAM: FIR LINEAR PHASE FILTER DESIGN PROGRAM
  5. C
  6. C AUTHORS: JAMES H. MCCLELLAN
  7. C          DEPARTMENT OF ELECTRICAL ENGINEERING AND COMPUTER SCIENCE
  8. C          MASSACHUSETTS INSTITUTE OF TECHNOLOGY
  9. C          CAMBRIDGE, MASS. 02139
  10. C
  11. C          THOMAS W. PARKS
  12. C          DEPARTMENT OF ELECTRICAL ENGINEERING
  13. C          RICE UNIVERSITY
  14. C          HOUSTON, TEXAS 77001
  15. C
  16. C          LAWRENCE R. RABINER
  17. C          BELL LABORATORIES
  18. C          MURRAY HILL, NEW JERSEY 07974
  19. C
  20. C INPUT:
  21. C  NFILT-- FILTER LENGTH
  22. C  JTYPE-- TYPE OF FILTER
  23. C          1 = MULTIPLE PASSBAND/STOPBAND FILTER
  24. C          2 = DIFFERENTIATOR
  25. C          3 = HILBERT TRANSFORM FILTER
  26. C  NBANDS-- NUMBER OF BANDS
  27. C  LGRID-- GRID DENSITY, WILL BE SET TO 16 UNLESS
  28. C          SPECIFIED OTHERWISE BY A POSITIVE CONSTANT.
  29. C
  30. C  EDGE(2*NBANDS)-- BANDEDGE ARRAY, LOWER AND UPPER EDGES FOR EACH BAND
  31. C                   WITH A MAXIMUM OF 10 BANDS.
  32. C
  33. C  FX(NBANDS)-- DESIRED FUNCTION ARRAY (OR DESIRED SLOPE IF A
  34. C               DIFFERENTIATOR) FOR EACH BAND.
  35. C
  36. C  WTX(NBANDS)-- WEIGHT FUNCTION ARRAY IN EACH BAND.  FOR A
  37. C                DIFFERENTIATOR, THE WEIGHT FUNCTION IS INVERSELY
  38. C                PROPORTIONAL TO F.
  39. C
  40. C  SAMPLE INPUT DATA SETUP:
  41. C       32,1,3,0
  42. C       0.0,0.1,0.2,0.35
  43. C       0.425,0.5
  44. C       0.0,1.0,0.0
  45. C       10.0,1.0,10.0
  46. C     THIS DATA SPECIFIES A LENGTH 32 BANDPASS FILTER WITH
  47. C     STOPBANDS 0 TO 0.1 AND 0.425 TO 0.5, AND PASSBAND FROM
  48. C     0.2 TO 0.35 WITH WEIGHTING OF 10 IN THE STOPBANDS AND 1
  49. C     IN THE PASSBAND.  THE GRID DENSITY DEFAULTS TO 16.
  50. C     THIS IS THE FILTER IN FIGURE 10.
  51. C
  52. C     THE FOLLOWING INPUT DATA SPECIFIES A LENGTH 32 FULLBAND
  53. C     DIFFERENTIATOR WITH SLOPE 1 AND WEIGHTING OF 1/F.
  54. C     THE GRID DENSITY WILL BE SET TO 20.
  55. C       32,2,1,20
  56. C       0,0.5
  57. C       1.0
  58. C       1.0
  59. C
  60. C-----------------------------------------------------------------------
  61. C
  62.       COMMON PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
  63.       COMMON /OOPS/NITER,IOUT
  64.       DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
  65.       DIMENSION H(66)
  66.       DIMENSION DES(1045),GRID(1045),WT(1045)
  67.       DIMENSION EDGE(20),FX(10),WTX(10),DEVIAT(10)
  68.       DOUBLE PRECISION PI2,PI
  69.       DOUBLE PRECISION AD,DEV,X,Y
  70.       DOUBLE PRECISION GEE,D
  71.       INTEGER BD1,BD2,BD3,BD4
  72.       DATA BD1,BD2,BD3,BD4/'B','A','N','D'/
  73. C     IOUT=1
  74. C     INPUT=2
  75.       PI=4.0*DATAN(1.0D0)
  76.       PI2=2.0D00*PI
  77. C
  78. C  THE PROGRAM IS SET UP FOR A MAXIMUM LENGTH OF 128, BUT
  79. C  THIS UPPER LIMIT CAN BE CHANGED BY REDIMENSIONING THE
  80. C  ARRAYS IEXT, AD, ALPHA, X, Y, H TO BE NFMAX/2 + 2.
  81. C  THE ARRAYS DES, GRID, AND WT MUST DIMENSIONED
  82. C  16(NFMAX/2 + 2).
  83. C
  84.       NFMAX=128
  85.   100 CONTINUE
  86.       JTYPE=0
  87. C
  88. C  PROGRAM INPUT SECTION
  89. C
  90.       READ(*,110) NFILT,JTYPE,NBANDS,LGRID
  91.       IF(NFILT.EQ.0)STOP
  92.   110 FORMAT(4I5)
  93.       IF(NFILT.LE.NFMAX.AND.NFILT.GE.3) GO TO 115
  94.       CALL ERROR
  95.       STOP
  96.   115 IF(NBANDS.LE.0) NBANDS=1
  97. C
  98. C  GRID DENSITY IS ASSUMED TO BE 16 UNLESS SPECIFIED
  99. C  OTHERWISE
  100. C
  101.       IF(LGRID.LE.0) LGRID=16
  102.       JB=2*NBANDS
  103.       READ(*,120) (EDGE(J),J=1,JB)
  104.   120 FORMAT(4F15.9)
  105.       READ(*,120) (FX(J),J=1,NBANDS)
  106.       READ(*,120) (WTX(J),J=1,NBANDS)
  107.       IF(JTYPE.GT.0.AND.JTYPE.LE.3) GO TO 125
  108.       CALL ERROR
  109.       STOP
  110.   125 NEG=1
  111.       IF(JTYPE.EQ.1) NEG=0
  112.       NODD=NFILT/2
  113.       NODD=NFILT-2*NODD
  114.       NFCNS=NFILT/2
  115.       IF(NODD.EQ.1.AND.NEG.EQ.0) NFCNS=NFCNS+1
  116. C
  117. C  SET UP THE DENSE GRID.  THE NUMBER OF POINTS IN THE GRID
  118. C  IS (FILTER LENGTH + 1)*GRID DENSITY/2
  119. C
  120.       GRID(1)=EDGE(1)
  121.       DELF=LGRID*NFCNS
  122.       DELF=0.5/DELF
  123.       IF(NEG.EQ.0) GO TO 135
  124.       IF(EDGE(1).LT.DELF) GRID(1)=DELF
  125.   135 CONTINUE
  126.       J=1
  127.       L=1
  128.       LBAND=1
  129.   140 FUP=EDGE(L+1)
  130.   145 TEMP=GRID(J)
  131. C
  132. C  CALCULATE THE DESIRED MAGNITUDE RESPONSE AND THE WEIGHT
  133. C  FUNCTION ON THE GRID
  134. C
  135.       DES(J)=EFF(TEMP,FX,WTX,LBAND,JTYPE)
  136.       WT(J)=WATE(TEMP,FX,WTX,LBAND,JTYPE)
  137.       J=J+1
  138.       GRID(J)=TEMP+DELF
  139.       IF(GRID(J).GT.FUP) GO TO 150
  140.       GO TO 145
  141.   150 GRID(J-1)=FUP
  142.       DES(J-1)=EFF(FUP,FX,WTX,LBAND,JTYPE)
  143.       WT(J-1)=WATE(FUP,FX,WTX,LBAND,JTYPE)
  144.       LBAND=LBAND+1
  145.       L=L+2
  146.       IF(LBAND.GT.NBANDS) GO TO 160
  147.       GRID(J)=EDGE(L)
  148.       GO TO 140
  149.   160 NGRID=J-1
  150.       IF(NEG.NE.NODD) GO TO 165
  151.       IF(GRID(NGRID).GT.(0.5-DELF)) NGRID=NGRID-1
  152.   165 CONTINUE
  153. C
  154. C  SET UP A NEW APPROXIMATION PROBLEM WHICH IS EQUIVALENT
  155. C  TO THE ORIGINAL PROBLEM
  156. C
  157.       IF(NEG) 170,170,180
  158.   170 IF(NODD.EQ.1) GO TO 200
  159.       DO 175 J=1,NGRID
  160.       CHANGE=DCOS(PI*GRID(J))
  161.       DES(J)=DES(J)/CHANGE
  162.   175 WT(J)=WT(J)*CHANGE
  163.       GO TO 200
  164.   180 IF(NODD.EQ.1) GO TO 190
  165.       DO 185 J=1,NGRID
  166.       CHANGE=DSIN(PI*GRID(J))
  167.       DES(J)=DES(J)/CHANGE
  168.   185 WT(J)=WT(J)*CHANGE
  169.       GO TO 200
  170.   190 DO 195 J=1,NGRID
  171.       CHANGE=DSIN(PI2*GRID(J))
  172.       DES(J)=DES(J)/CHANGE
  173.   195 WT(J)=WT(J)*CHANGE
  174. C
  175. C  INITIAL GUESS FOR THE EXTREMAL FREQUENCIES--EQUALLY
  176. C  SPACED ALONG THE GRID
  177. C
  178.   200 TEMP=FLOAT(NGRID-1)/FLOAT(NFCNS)
  179.       DO 210 J=1,NFCNS
  180.       XT=J-1
  181.   210 IEXT(J)=XT*TEMP+1.0
  182.       IEXT(NFCNS+1)=NGRID
  183.       NM1=NFCNS-1
  184.       NZ=NFCNS+1
  185. C
  186. C  CALL THE REMEZ EXCHANGE ALGORITHM TO DO THE APPROXIMATION
  187. C  PROBLEM
  188. C
  189.       CALL REMEZ
  190. C
  191. C  CALCULATE THE IMPULSE RESPONSE.
  192. C
  193.       IF(NEG) 300,300,320
  194.   300 IF(NODD.EQ.0) GO TO 310
  195.       DO 305 J=1,NM1
  196.       NZMJ=NZ-J
  197.   305 H(J)=0.5*ALPHA(NZMJ)
  198.       H(NFCNS)=ALPHA(1)
  199.       GO TO 350
  200.   310 H(1)=0.25*ALPHA(NFCNS)
  201.       DO 315 J=2,NM1
  202.       NZMJ=NZ-J
  203.       NF2J=NFCNS+2-J
  204.   315 H(J)=0.25*(ALPHA(NZMJ)+ALPHA(NF2J))
  205.       H(NFCNS)=0.5*ALPHA(1)+0.25*ALPHA(2)
  206.       GO TO 350
  207.   320 IF(NODD.EQ.0) GO TO 330
  208.       H(1)=0.25*ALPHA(NFCNS)
  209.       H(2)=0.25*ALPHA(NM1)
  210.       DO 325 J=3,NM1
  211.       NZMJ=NZ-J
  212.       NF3J=NFCNS+3-J
  213.   325 H(J)=0.25*(ALPHA(NZMJ)-ALPHA(NF3J))
  214.       H(NFCNS)=0.5*ALPHA(1)-0.25*ALPHA(3)
  215.       H(NZ)=0.0
  216.       GO TO 350
  217.   330 H(1)=0.25*ALPHA(NFCNS)
  218.       DO 335 J=2,NM1
  219.       NZMJ=NZ-J
  220.       NF2J=NFCNS+2-J
  221.   335 H(J)=0.25*(ALPHA(NZMJ)-ALPHA(NF2J))
  222.       H(NFCNS)=0.5*ALPHA(1)-0.25*ALPHA(2)
  223. C
  224. C  PROGRAM OUTPUT SECTION.
  225. C
  226.   350 WRITE(*,360)
  227.   360 FORMAT(1H , 70(1H*)//15X,29HFINITE IMPULSE RESPONSE (FIR)/
  228.      113X,34HLINEAR PHASE DIGITAL FILTER DESIGN/
  229.      217X,24HREMEZ EXCHANGE ALGORITHM/)
  230.       IF(JTYPE.EQ.1) WRITE(*,365)
  231.   365 FORMAT(22X,15HBANDPASS FILTER/)
  232.       IF(JTYPE.EQ.2) WRITE(*,370)
  233.   370 FORMAT(22X,14HDIFFERENTIATOR/)
  234.       IF(JTYPE.EQ.3) WRITE(*,375)
  235.   375 FORMAT(20X,19HHILBERT TRANSFORMER/)
  236.       WRITE(*,378) NFILT
  237.   378 FORMAT(20X,16HFILTER LENGTH = ,I3/)
  238.       WRITE(*,380)
  239.   380 FORMAT(15X,28H***** IMPULSE RESPONSE *****)
  240.       DO 381 J=1,NFCNS
  241.       K=NFILT+1-J
  242.       XFUK=H(j)*32767.0
  243.       IF(NEG.EQ.0) WRITE(*,382) J,xfuk,K
  244.       IF(NEG.EQ.1) WRITE(*,383) J,xfuk,K
  245.   381 CONTINUE
  246.   382 FORMAT(10X,2HH(,I2,4H) = ,F11.3,5H = H(,I3,1H))
  247.   383 FORMAT(10X,2HH(,I2,4H) = ,F11.3,6H = -H(,I3,1H))
  248.       IF(NEG.EQ.1.AND.NODD.EQ.1) WRITE(*,384) NZ
  249.   384 FORMAT(13X,2HH(,I2,8H) =  0.0)
  250.       DO 450 K=1,NBANDS,4
  251.       KUP=K+3
  252.       IF(KUP.GT.NBANDS) KUP=NBANDS
  253.       WRITE(*,385) (BD1,BD2,BD3,BD4,J,J=K,KUP)
  254.   385 FORMAT(/24X,4(4A1,I3,7X))
  255.       WRITE(*,390) (EDGE(2*J-1),J=K,KUP)
  256.   390 FORMAT(2X,15HLOWER BAND EDGE,5F14.7)
  257.       WRITE(*,395) (EDGE(2*J),J=K,KUP)
  258.   395 FORMAT(2X,15HUPPER BAND EDGE,5F14.7)
  259.       IF(JTYPE.NE.2) WRITE(*,400) (FX(J),J=K,KUP)
  260.   400 FORMAT(2X,13HDESIRED VALUE,2X,5F14.7)
  261.       IF(JTYPE.EQ.2) WRITE(*,405) (FX(J),J=K,KUP)
  262.   405 FORMAT(2X,13HDESIRED SLOPE,2X,5F14.7)
  263.       WRITE(*,410) (WTX(J),J=K,KUP)
  264.   410 FORMAT(2X,9HWEIGHTING,6X,5F14.7)
  265.       DO 420 J=K,KUP
  266.   420 DEVIAT(J)=DEV/WTX(J)
  267.       WRITE(*,425) (DEVIAT(J),J=K,KUP)
  268.   425 FORMAT(2X,9HDEVIATION,6X,5F14.7)
  269.       IF(JTYPE.NE.1) GO TO 450
  270.       DO 430 J=K,KUP
  271.   430 DEVIAT(J)=20.0*ALOG10(DEVIAT(J)+FX(J))
  272.       WRITE(*,435) (DEVIAT(J),J=K,KUP)
  273.   435 FORMAT(2X,15HDEVIATION IN DB,5F14.7)
  274.   450 CONTINUE
  275.       DO 452 J=1,NZ
  276.       IX=IEXT(J)
  277.   452 GRID(J)=GRID(IX)
  278.       WRITE(*,455) (GRID(J),J=1,NZ)
  279.   455 FORMAT(/2X,47HEXTREMAL FREQUENCIES--MAXIMA OF THE ERROR CURVE/
  280.      1 (2X,5F12.7))
  281.       WRITE(*,460)
  282.   460 FORMAT(/1X,70(1H*))
  283.       GO TO 100
  284.       END
  285. C
  286. C-----------------------------------------------------------------------
  287. C FUNCTION: EFF
  288. C   FUNCTION TO CALCULATE THE DESIRED MAGNITUDE RESPONSE
  289. C   AS A FUNCTION OF FREQUENCY.
  290. C   AN ARBITRARY FUNCTION OF FREQUENCY CAN BE
  291. C   APPROXIMATED IF THE USER REPLACES THIS FUNCTION
  292. C   WITH THE APPROPRIATE CODE TO EVALUATE THE IDEAL
  293. C   MAGNITUDE.  NOTE THAT THE PARAMETER FREQ IS THE
  294. C   VALUE OF NORMALIZED FREQUENCY NEEDED FOR EVALUATION.
  295. C-----------------------------------------------------------------------
  296. C
  297.       FUNCTION EFF(FREQ,FX,WTX,LBAND,JTYPE)
  298.       DIMENSION FX(5),WTX(5)
  299.       IF(JTYPE.EQ.2) GO TO 1
  300.       EFF=FX(LBAND)
  301.       RETURN
  302.     1 EFF=FX(LBAND)*FREQ
  303.       RETURN
  304.       END
  305. C
  306. C-----------------------------------------------------------------------
  307. C FUNCTION: WATE
  308. C   FUNCTION TO CALCULATE THE WEIGHT FUNCTION AS A FUNCTION
  309. C   OF FREQUENCY.  SIMILAR TO THE FUNCTION EFF, THIS FUNCTION CAN
  310. C   BE REPLACED BY A USER-WRITTEN ROUTINE TO CALCULATE ANY
  311. C   DESIRED WEIGHTING FUNCTION.
  312. C-----------------------------------------------------------------------
  313. C
  314.       FUNCTION WATE(FREQ,FX,WTX,LBAND,JTYPE)
  315.       DIMENSION FX(5),WTX(5)
  316.       IF(JTYPE.EQ.2) GO TO 1
  317.       WATE=WTX(LBAND)
  318.       RETURN
  319.     1 IF(FX(LBAND).LT.0.0001) GO TO 2
  320.       WATE=WTX(LBAND)/FREQ
  321.       RETURN
  322.     2 WATE=WTX(LBAND)
  323.       RETURN
  324.       END
  325. C
  326. C-----------------------------------------------------------------------
  327. C SUBROUTINE: ERROR
  328. C   THIS ROUTINE WRITES AN ERROR MESSAGE IF AN
  329. C   ERROR HAS BEEN DETECTED IN THE INPUT DATA.
  330. C-----------------------------------------------------------------------
  331. C
  332.       SUBROUTINE ERROR
  333.       COMMON /OOPS/NITER,IOUT
  334.       WRITE(*,1)
  335.     1 FORMAT(44H ************ ERROR IN INPUT DATA **********)
  336.       RETURN
  337.       END
  338. C
  339. C-----------------------------------------------------------------------
  340. C SUBROUTINE: REMEZ
  341. C   THIS SUBROUTINE IMPLEMENTS THE REMEZ EXCHANGE ALGORITHM
  342. C   FOR THE WEIGHTED CHEBYSHEV APPROXIMATION OF A CONTINUOUS
  343. C   FUNCTION WITH A SUM OF COSINES.  INPUTS TO THE SUBROUTINE
  344. C   ARE A DENSE GRID WHICH REPLACES THE FREQUENCY AXIS, THE
  345. C   DESIRED FUNCTION ON THIS GRID, THE WEIGHT FUNCTION ON THE
  346. C   GRID, THE NUMBER OF COSINES, AND AN INITIAL GUESS OF THE
  347. C   EXTREMAL FREQUENCIES.  THE PROGRAM MINIMIZES THE CHEBYSHEV
  348. C   ERROR BY DETERMINING THE BEST LOCATION OF THE EXTREMAL
  349. C   FREQUENCIES (POINTS OF MAXIMUM ERROR) AND THEN CALCULATES
  350. C   THE COEFFICIENTS OF THE BEST APPROXIMATION.
  351. C-----------------------------------------------------------------------
  352. C
  353.       SUBROUTINE REMEZ
  354.       COMMON PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
  355.       COMMON /OOPS/NITER,IOUT
  356.       DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
  357.       DIMENSION DES(1045),GRID(1045),WT(1045)
  358.       DIMENSION A(66),P(65),Q(65)
  359.       DOUBLE PRECISION PI2,DNUM,DDEN,DTEMP,A,P,Q
  360.       DOUBLE PRECISION DK,DAK
  361.       DOUBLE PRECISION AD,DEV,X,Y
  362.       DOUBLE PRECISION GEE,D
  363. C
  364. C  THE PROGRAM ALLOWS A MAXIMUM NUMBER OF ITERATIONS OF 25
  365. C
  366.       ITRMAX=25
  367.       DEVL=-1.0
  368.       NZ=NFCNS+1
  369.       NZZ=NFCNS+2
  370.       NITER=0
  371.   100 CONTINUE
  372.       IEXT(NZZ)=NGRID+1
  373.       NITER=NITER+1
  374.       IF(NITER.GT.ITRMAX) GO TO 400
  375.       DO 110 J=1,NZ
  376.       JXT=IEXT(J)
  377.       DTEMP=GRID(JXT)
  378.       DTEMP=DCOS(DTEMP*PI2)
  379.   110 X(J)=DTEMP
  380.       JET=(NFCNS-1)/15+1
  381.       DO 120 J=1,NZ
  382.   120 AD(J)=D(J,NZ,JET)
  383.       DNUM=0.0
  384.       DDEN=0.0
  385.       K=1
  386.       DO 130 J=1,NZ
  387.       L=IEXT(J)
  388.       DTEMP=AD(J)*DES(L)
  389.       DNUM=DNUM+DTEMP
  390.       DTEMP=FLOAT(K)*AD(J)/WT(L)
  391.       DDEN=DDEN+DTEMP
  392.   130 K=-K
  393.       DEV=DNUM/DDEN
  394.       WRITE(*,131) DEV
  395.   131 FORMAT(1X,12HDEVIATION = ,F12.9)
  396.       NU=1
  397.       IF(DEV.GT.0.0) NU=-1
  398.       DEV=-FLOAT(NU)*DEV
  399.       K=NU
  400.       DO 140 J=1,NZ
  401.       L=IEXT(J)
  402.       DTEMP=FLOAT(K)*DEV/WT(L)
  403.       Y(J)=DES(L)+DTEMP
  404.   140 K=-K
  405.       IF(DEV.GT.DEVL) GO TO 150
  406.       CALL OUCH
  407.       GO TO 400
  408.   150 DEVL=DEV
  409.       JCHNGE=0
  410.       K1=IEXT(1)
  411.       KNZ=IEXT(NZ)
  412.       KLOW=0
  413.       NUT=-NU
  414.       J=1
  415. C
  416. C  SEARCH FOR THE EXTREMAL FREQUENCIES OF THE BEST
  417. C  APPROXIMATION
  418. C
  419.   200 IF(J.EQ.NZZ) YNZ=COMP
  420.       IF(J.GE.NZZ) GO TO 300
  421.       KUP=IEXT(J+1)
  422.       L=IEXT(J)+1
  423.       NUT=-NUT
  424.       IF(J.EQ.2) Y1=COMP
  425.       COMP=DEV
  426.       IF(L.GE.KUP) GO TO 220
  427.       ERR=GEE(L,NZ)
  428.       ERR=(ERR-DES(L))*WT(L)
  429.       DTEMP=FLOAT(NUT)*ERR-COMP
  430.       IF(DTEMP.LE.0.0) GO TO 220
  431.       COMP=FLOAT(NUT)*ERR
  432.   210 L=L+1
  433.       IF(L.GE.KUP) GO TO 215
  434.       ERR=GEE(L,NZ)
  435.       ERR=(ERR-DES(L))*WT(L)
  436.       DTEMP=FLOAT(NUT)*ERR-COMP
  437.       IF(DTEMP.LE.0.0) GO TO 215
  438.       COMP=FLOAT(NUT)*ERR
  439.       GO TO 210
  440.   215 IEXT(J)=L-1
  441.       J=J+1
  442.       KLOW=L-1
  443.       JCHNGE=JCHNGE+1
  444.       GO TO 200
  445.   220 L=L-1
  446.   225 L=L-1
  447.       IF(L.LE.KLOW) GO TO 250
  448.       ERR=GEE(L,NZ)
  449.       ERR=(ERR-DES(L))*WT(L)
  450.       DTEMP=FLOAT(NUT)*ERR-COMP
  451.       IF(DTEMP.GT.0.0) GO TO 230
  452.       IF(JCHNGE.LE.0) GO TO 225
  453.       GO TO 260
  454.   230 COMP=FLOAT(NUT)*ERR
  455.   235 L=L-1
  456.       IF(L.LE.KLOW) GO TO 240
  457.       ERR=GEE(L,NZ)
  458.       ERR=(ERR-DES(L))*WT(L)
  459.       DTEMP=FLOAT(NUT)*ERR-COMP
  460.       IF(DTEMP.LE.0.0) GO TO 240
  461.       COMP=FLOAT(NUT)*ERR
  462.       GO TO 235
  463.   240 KLOW=IEXT(J)
  464.       IEXT(J)=L+1
  465.       J=J+1
  466.       JCHNGE=JCHNGE+1
  467.       GO TO 200
  468.   250 L=IEXT(J)+1
  469.       IF(JCHNGE.GT.0) GO TO 215
  470.   255 L=L+1
  471.       IF(L.GE.KUP) GO TO 260
  472.       ERR=GEE(L,NZ)
  473.       ERR=(ERR-DES(L))*WT(L)
  474.       DTEMP=FLOAT(NUT)*ERR-COMP
  475.       IF(DTEMP.LE.0.0) GO TO 255
  476.       COMP=FLOAT(NUT)*ERR
  477.       GO TO 210
  478.   260 KLOW=IEXT(J)
  479.       J=J+1
  480.       GO TO 200
  481.   300 IF(J.GT.NZZ) GO TO 320
  482.       IF(K1.GT.IEXT(1)) K1=IEXT(1)
  483.       IF(KNZ.LT.IEXT(NZ)) KNZ=IEXT(NZ)
  484.       NUT1=NUT
  485.       NUT=-NU
  486.       L=0
  487.       KUP=K1
  488.       COMP=YNZ*(1.00001)
  489.       LUCK=1
  490.   310 L=L+1
  491.       IF(L.GE.KUP) GO TO 315
  492.       ERR=GEE(L,NZ)
  493.       ERR=(ERR-DES(L))*WT(L)
  494.       DTEMP=FLOAT(NUT)*ERR-COMP
  495.       IF(DTEMP.LE.0.0) GO TO 310
  496.       COMP=FLOAT(NUT)*ERR
  497.       J=NZZ
  498.       GO TO 210
  499.   315 LUCK=6
  500.       GO TO 325
  501.   320 IF(LUCK.GT.9) GO TO 350
  502.       IF(COMP.GT.Y1) Y1=COMP
  503.       K1=IEXT(NZZ)
  504.   325 L=NGRID+1
  505.       KLOW=KNZ
  506.       NUT=-NUT1
  507.       COMP=Y1*(1.00001)
  508.   330 L=L-1
  509.       IF(L.LE.KLOW) GO TO 340
  510.       ERR=GEE(L,NZ)
  511.       ERR=(ERR-DES(L))*WT(L)
  512.       DTEMP=FLOAT(NUT)*ERR-COMP
  513.       IF(DTEMP.LE.0.0) GO TO 330
  514.       J=NZZ
  515.       COMP=FLOAT(NUT)*ERR
  516.       LUCK=LUCK+10
  517.       GO TO 235
  518.   340 IF(LUCK.EQ.6) GO TO 370
  519.       DO 345 J=1,NFCNS
  520.       NZZMJ=NZZ-J
  521.       NZMJ=NZ-J
  522.   345 IEXT(NZZMJ)=IEXT(NZMJ)
  523.       IEXT(1)=K1
  524.       GO TO 100
  525.   350 KN=IEXT(NZZ)
  526.       DO 360 J=1,NFCNS
  527.   360 IEXT(J)=IEXT(J+1)
  528.       IEXT(NZ)=KN
  529.       GO TO 100
  530.   370 IF(JCHNGE.GT.0) GO TO 100
  531. C
  532. C  CALCULATION OF THE COEFFICIENTS OF THE BEST APPROXIMATION
  533. C  USING THE INVERSE DISCRETE FOURIER TRANSFORM
  534. C
  535.   400 CONTINUE
  536.       NM1=NFCNS-1
  537.       FSH=1.0E-06
  538.       GTEMP=GRID(1)
  539.       X(NZZ)=-2.0
  540.       CN=2*NFCNS-1
  541.       DELF=1.0/CN
  542.       L=1
  543.       KKK=0
  544.       IF(GRID(1).LT.0.01.AND.GRID(NGRID).GT.0.49) KKK=1
  545.       IF(NFCNS.LE.3) KKK=1
  546.       IF(KKK.EQ.1) GO TO 405
  547.       DTEMP=DCOS(PI2*GRID(1))
  548.       DNUM=DCOS(PI2*GRID(NGRID))
  549.       AA=2.0/(DTEMP-DNUM)
  550.       BB=-(DTEMP+DNUM)/(DTEMP-DNUM)
  551.   405 CONTINUE
  552.       DO 430 J=1,NFCNS
  553.       FT=J-1
  554.       FT=FT*DELF
  555.       XT=DCOS(PI2*FT)
  556.       IF(KKK.EQ.1) GO TO 410
  557.       XT=(XT-BB)/AA
  558.       XT1=SQRT(1.0-XT*XT)
  559.       FT=ATAN2(XT1,XT)/PI2
  560.   410 XE=X(L)
  561.       IF(XT.GT.XE) GO TO 420
  562.       IF((XE-XT).LT.FSH) GO TO 415
  563.       L=L+1
  564.       GO TO 410
  565.   415 A(J)=Y(L)
  566.       GO TO 425
  567.   420 IF((XT-XE).LT.FSH) GO TO 415
  568.       GRID(1)=FT
  569.       A(J)=GEE(1,NZ)
  570.   425 CONTINUE
  571.       IF(L.GT.1) L=L-1
  572.   430 CONTINUE
  573.       GRID(1)=GTEMP
  574.       DDEN=PI2/CN
  575.       DO 510 J=1,NFCNS
  576.       DTEMP=0.0
  577.       DNUM=J-1
  578.       DNUM=DNUM*DDEN
  579.       IF(NM1.LT.1) GO TO 505
  580.       DO 500 K=1,NM1
  581.       DAK=A(K+1)
  582.       DK=K
  583.   500 DTEMP=DTEMP+DAK*DCOS(DNUM*DK)
  584.   505 DTEMP=2.0*DTEMP+A(1)
  585.   510 ALPHA(J)=DTEMP
  586.       DO 550 J=2,NFCNS
  587.   550 ALPHA(J)=2.0*ALPHA(J)/CN
  588.       ALPHA(1)=ALPHA(1)/CN
  589.       IF(KKK.EQ.1) GO TO 545
  590.       P(1)=2.0*ALPHA(NFCNS)*BB+ALPHA(NM1)
  591.       P(2)=2.0*AA*ALPHA(NFCNS)
  592.       Q(1)=ALPHA(NFCNS-2)-ALPHA(NFCNS)
  593.       DO 540 J=2,NM1
  594.       IF(J.LT.NM1) GO TO 515
  595.       AA=0.5*AA
  596.       BB=0.5*BB
  597.   515 CONTINUE
  598.       P(J+1)=0.0
  599.       DO 520 K=1,J
  600.       A(K)=P(K)
  601.   520 P(K)=2.0*BB*A(K)
  602.       P(2)=P(2)+A(1)*2.0*AA
  603.       JM1=J-1
  604.       DO 525 K=1,JM1
  605.   525 P(K)=P(K)+Q(K)+AA*A(K+1)
  606.       JP1=J+1
  607.       DO 530 K=3,JP1
  608.   530 P(K)=P(K)+AA*A(K-1)
  609.       IF(J.EQ.NM1) GO TO 540
  610.       DO 535 K=1,J
  611.   535 Q(K)=-A(K)
  612.       NF1J=NFCNS-1-J
  613.       Q(1)=Q(1)+ALPHA(NF1J)
  614.   540 CONTINUE
  615.       DO 543 J=1,NFCNS
  616.   543 ALPHA(J)=P(J)
  617.   545 CONTINUE
  618.       IF(NFCNS.GT.3) RETURN
  619.       ALPHA(NFCNS+1)=0.0
  620.       ALPHA(NFCNS+2)=0.0
  621.       RETURN
  622.       END
  623. C
  624. C-----------------------------------------------------------------------
  625. C FUNCTION: D
  626. C   FUNCTION TO CALCULATE THE LAGRANGE INTERPOLATION
  627. C   COEFFICIENTS FOR USE IN THE FUNCTION GEE.
  628. C-----------------------------------------------------------------------
  629. C
  630.       DOUBLE PRECISION FUNCTION D(K,N,M)
  631.       COMMON PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
  632.       DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
  633.       DIMENSION DES(1045),GRID(1045),WT(1045)
  634.       DOUBLE PRECISION AD,DEV,X,Y
  635.       DOUBLE PRECISION Q
  636.       DOUBLE PRECISION PI2
  637.       D=1.0
  638.       Q=X(K)
  639.       DO 3 L=1,M
  640.       DO 2 J=L,N,M
  641.       IF(J-K)1,2,1
  642.     1 D=2.0*D*(Q-X(J))
  643.     2 CONTINUE
  644.     3 CONTINUE
  645.       D=1.0/D
  646.       RETURN
  647.       END
  648. C
  649. C-----------------------------------------------------------------------
  650. C FUNCTION: GEE
  651. C   FUNCTION TO EVALUATE THE FREQUENCY RESPONSE USING THE
  652. C   LAGRANGE INTERPOLATION FORMULA IN THE BARYCENTRIC FORM
  653. C-----------------------------------------------------------------------
  654. C
  655.       DOUBLE PRECISION FUNCTION GEE(K,N)
  656.       COMMON PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
  657.       DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
  658.       DIMENSION DES(1045),GRID(1045),WT(1045)
  659.       DOUBLE PRECISION P,C,D,XF
  660.       DOUBLE PRECISION PI2
  661.       DOUBLE PRECISION AD,DEV,X,Y
  662.       P=0.0
  663.       XF=GRID(K)
  664.       XF=DCOS(PI2*XF)
  665.       D=0.0
  666.       DO 1 J=1,N
  667.       C=XF-X(J)
  668.       C=AD(J)/C
  669.       D=D+C
  670.     1 P=P+C*Y(J)
  671.       GEE=P/D
  672.       RETURN
  673.       END
  674. C
  675. C-----------------------------------------------------------------------
  676. C SUBROUTINE: OUCH
  677. C   WRITES AN ERROR MESSAGE WHEN THE ALGORITHM FAILS TO
  678. C   CONVERGE.  THERE SEEM TO BE TWO CONDITIONS UNDER WHICH
  679. C   THE ALGORITHM FAILS TO CONVERGE: (1) THE INITIAL
  680. C   GUESS FOR THE EXTREMAL FREQUENCIES IS SO POOR THAT
  681. C   THE EXCHANGE ITERATION CANNOT GET STARTED, OR
  682. C   (2) NEAR THE TERMINATION OF A CORRECT DESIGN,
  683. C   THE DEVIATION DECREASES DUE TO ROUNDING ERRORS
  684. C   AND THE PROGRAM STOPS.  IN THIS LATTER CASE THE
  685. C   FILTER DESIGN IS PROBABLY ACCEPTABLE, BUT SHOULD
  686. C   BE CHECKED BY COMPUTING A FREQUENCY RESPONSE.
  687. C-----------------------------------------------------------------------
  688. C
  689.       SUBROUTINE OUCH
  690.       COMMON /OOPS/NITER,IOUT
  691.       WRITE(*,1)NITER
  692.     1 FORMAT(44H ************ FAILURE TO CONVERGE **********/
  693.      141H0PROBABLE CAUSE IS MACHINE ROUNDING ERROR/
  694.      223H0NUMBER OF ITERATIONS =,I4/
  695.      339H0IF THE NUMBER OF ITERATIONS EXCEEDS 3,/
  696.      462H0THE DESIGN MAY BE CORRECT, BUT SHOULD BE VERIFIED WITH AN FFT)
  697.       RETURN
  698.       END
  699.  
  700.  
  701.